#
# CCES2016.r
#
rm(list=ls(all=TRUE))
#
library(haven)
library(plyr)
library(dplyr)
library(questionr)
library(statar)
library(ggplot2)
library(viridis)
library(rjags)
library(apsrtable)
library(mokken)
library(basicspace)
library(gridExtra)
library(psych)
library(interflex)
library(ggpubr)
#
setwd("c:/")
data <- read_dta("Dropbox/Networks/data/CCES16_UCDUGA_withcounty.dta")
#
attach(data,warn.conflicts = FALSE)
#
rescale01 <- function(x){(x-min(x, na.rm=TRUE))/(max(x, na.rm=TRUE)-min(x, na.rm=TRUE))}
higher.conservative <- function(x){
	if (cor(x, ideology, use="complete") < 0) {
		(x <- -1 * x + max(x, na.rm=TRUE))}
	else{x <- x}}
#
#  PARTY ID (1 = Democrat, 2 = Independent, 3 = Republican)
#
partyid <- NA
partyid[pid7==1 | pid7==2 | pid7==3] <- 1
partyid[pid7==4] <- 2
partyid[pid7==5 | pid7==6 | pid7==7] <- 3
#
partyid7 <- pid7
partyid7[partyid7 < 1 | partyid7 > 7] <- NA
#
partisanextremity <- abs(partyid7 - 4)
#
#
# PRESIDENTIAL VOTE
#
presvote <- NA
presvote[CC16_410a==2] <- "Clinton"
presvote[CC16_410a==1] <- "Trump"
presvote[CC16_410a==3 | CC16_410a==8] <- "Other/not sure"
presvote <- as.factor(presvote)
#
#
# DEMOGRAPHICS
#
education <- NA
education[educ==1 | educ==2] <- 1
education[educ==3 | educ==4] <- 2
education[educ==5 | educ==6] <- 3
education <- rescale01(education)
#
income <- NA
income[faminc>=1 & faminc <=4] <- 1
income[faminc>=5 & faminc <=8] <- 2
income[faminc>=9 & faminc <=31] <- 3
income <- rescale01(income)
#
age <- 2016 - birthyr
age[age>=18 & age <=30] <- 1
age[age>=31 & age <=45] <- 2
age[age>=46 & age <=60] <- 3
age[age>=61] <- 4
age <- rescale01(age)
#
female <- gender - 1
female[female < 0 | female > 1] <- NA
#
black <- NA
black[race!=2] <- 1
black[race==2] <- 2
#
latino <- NA
latino[race==3 | hispanic==1] <- 2
latino[race!=3 & hispanic!=1] <- 1
#
unionmember <- NA
unionmember[union==1 | union==2 | unionhh==1 | unionhh==2] <- 2
unionmember[union==3 & unionhh==3] <- 1
unionmember[unionmember > 2] <- NA
#
bornagain <- pew_bornagain
bornagain <- -1 * (bornagain - 3)
bornagain[bornagain > 2] <- NA
bornagain <- (bornagain-min(bornagain, na.rm=TRUE))/(max(bornagain, na.rm=TRUE)-min(bornagain, na.rm=TRUE))
#
importancereligion <- pew_religimp
importancereligion <- -1 * (importancereligion - 5)
importancereligion[importancereligion > 4] <- NA
importancereligion <- (importancereligion-min(importancereligion, na.rm=TRUE))/(max(importancereligion, na.rm=TRUE)-min(importancereligion, na.rm=TRUE))
#
churchattendance <- pew_churatd
churchattendance[churchattendance > 6] <- NA
churchattendance <- -1 * (churchattendance - 7)
churchattendance <- (churchattendance-min(churchattendance, na.rm=TRUE))/(max(churchattendance, na.rm=TRUE)-min(churchattendance, na.rm=TRUE))
#
frequencyprayer <- pew_prayer
frequencyprayer[frequencyprayer > 7] <- NA
frequencyprayer <- -1 * (frequencyprayer - 8)
frequencyprayer <- (frequencyprayer-min(frequencyprayer, na.rm=TRUE))/(max(frequencyprayer, na.rm=TRUE)-min(frequencyprayer, na.rm=TRUE))
#
religiositymat <- cbind(bornagain, importancereligion, churchattendance, frequencyprayer)
psych::alpha(religiositymat)
religiosity.scores <- rowMeans(religiositymat, na.rm=TRUE)
religiosity.binary <- ntile(religiosity.scores, 3)
#
rr1 <- CC16_422c
rr1[rr1 < 1 | rr1 > 5] <- NA
#
rr2 <- CC16_422d
rr2[rr2 < 1 | rr2 > 5] <- NA
#
rr3 <- CC16_422e
rr3[rr3 < 1 | rr3 > 5] <- NA
rr3 <- -1 * (rr3 - 6)
#
rr4 <- CC16_422f
rr4[rr4 < 1 | rr4 > 5] <- NA
rr4 <- -1 * (rr4 - 6)
#
rrmat <- cbind(rr1, rr2, rr4)
psych::alpha(rrmat)
racialresentment.scores <- rowMeans(rrmat, na.rm=TRUE)
racialresentment.binary <- ntile(racialresentment.scores, 3)
#
#
# COUNTY-LEVEL INFORMATION
#
county.ideology <- mrp_ideology_mean
#
county.gop2016 <- per_gop_2016
#
county.delta.gop2012gop2016 <- per_gop_2016 - per_gop_2012
#
county.trumpprimary2016 <- trump_2016primary
#
county.sandersprimary2016 <- sanders_2016primary
#
#
# IDEOLOGICAL SELF-PLACEMENT
#
ideological.selfplace <- CC16_340a
ideological.selfplace[ideological.selfplace < 1 | ideological.selfplace > 7] <- NA
#
#
# SUMMATED IDEOLOGICAL SCALE
#
ideology <- CC16_340a
ideology[ideology < 1 | ideology > 7] <- NA
#
raiseminimumwage <- CC16_351K
raiseminimumwage[raiseminimumwage > 2] <- NA
#
repealACA <- CC16_351I
repealACA[repealACA > 2] <- NA
#
regulateCO2 <- CC16_333a
regulateCO2[regulateCO2 > 2] <- NA
#
renewablefuels <- CC16_333c
renewablefuels[renewablefuels > 2] <- NA
#
enforcementcleanairact <- CC16_333d
enforcementcleanairact[enforcementcleanairact > 2] <- NA
#
syrianrefugees <- CC16_331_5
syrianrefugees[syrianrefugees > 2] <- NA
#
muslimban <- CC16_331_8
muslimban[muslimban > 2] <- NA
#
isis_sendsignificantforce <- CC16_312_7
isis_sendsignificantforce[isis_sendsignificantforce > 2] <- NA
#
immigration_grantlegalstatus <- CC16_331_1
immigration_grantlegalstatus[immigration_grantlegalstatus > 2] <- NA
#
immigration_increasepatrols <- CC16_331_2
immigration_increasepatrols <- -1 * (immigration_increasepatrols - 3)
immigration_increasepatrols[immigration_increasepatrols > 2] <- NA
#
immigration_dreamers <- CC16_331_3
immigration_dreamers[immigration_dreamers > 2] <- NA
#
immigration_finebusinesseshire <- CC16_331_4
immigration_finebusinesseshire <- -1 * (immigration_finebusinesseshire - 3)
immigration_finebusinesseshire[immigration_finebusinesseshire > 2] <- NA
#
immigration_deportillegalimmigrants <- CC16_331_7
immigration_deportillegalimmigrants <- -1 * (immigration_deportillegalimmigrants - 3)
immigration_deportillegalimmigrants[immigration_deportillegalimmigrants > 2] <- NA
#
guncontrol_banassaultrifles <- CC16_330d
guncontrol_banassaultrifles[guncontrol_banassaultrifles > 2] <- NA
#
guncontrol_easierconcealcarry <- CC16_330e
guncontrol_easierconcealcarry[guncontrol_easierconcealcarry > 2] <- NA
#
alwaysallowabortion <- CC16_332a
alwaysallowabortion[alwaysallowabortion > 2] <- NA
#
prohibitabortiontwentyweeks <- CC16_332c
prohibitabortiontwentyweeks[prohibitabortiontwentyweeks > 2] <- NA
prohibitabortiontwentyweeks <- (prohibitabortiontwentyweeks - 3) * -1
#
allowemployersdeclineabortion <- CC16_332d
allowemployersdeclineabortion[allowemployersdeclineabortion > 2] <- NA
allowemployersdeclineabortion <- (allowemployersdeclineabortion - 3) * -1
#
prohibitfederalfundingabortions <- CC16_332e
prohibitfederalfundingabortions[prohibitfederalfundingabortions > 2] <- NA
prohibitfederalfundingabortions <- (prohibitfederalfundingabortions - 3) * -1
#
gaymarriage <- CC16_335
gaymarriage[gaymarriage > 2] <- NA
#
issues <- cbind(raiseminimumwage, repealACA, regulateCO2,
	renewablefuels, enforcementcleanairact, syrianrefugees,
	muslimban, isis_sendsignificantforce, immigration_grantlegalstatus,
	immigration_increasepatrols,
	immigration_dreamers, immigration_finebusinesseshire,
	immigration_deportillegalimmigrants,
	guncontrol_banassaultrifles,
	guncontrol_easierconcealcarry, alwaysallowabortion, prohibitabortiontwentyweeks,
	allowemployersdeclineabortion, prohibitfederalfundingabortions, gaymarriage)
#
issues <- apply(issues, 2, rescale01)
issues <- apply(issues, 2, higher.conservative)
#
scale <- aisp(na.omit(issues), lowerbound=0.5)
#
policypreferences <- scale(rowMeans(issues[,scale==1], na.rm=TRUE))
#
#
#
#
# NUMBER OF DISCUSSANTS
#
numberdiscussants <-  NA
numberdiscussants[legitimate_name_1==0] <- 0
numberdiscussants[!is.na(UCD401_1) & legitimate_name_1==1] <- 1
numberdiscussants[!is.na(UCD401_2) & legitimate_name_1==1] <- 2
numberdiscussants[!is.na(UCD401_3) & legitimate_name_1==1] <- 3
numberdiscussants[is.na(numberdiscussants)] <- 0
#
#
#  DISCUSSANT VOTE (1 = Clinton, 2 = Trump, 3 = Other/Not Sure, 0 = Missing)
#
discussant.1.vote <- NA
discussant.1.vote[UCD405A==1 & numberdiscussants >= 1] <- 1
discussant.1.vote[UCD405A==2 & numberdiscussants >= 1] <- 2
discussant.1.vote[(UCD405A==3 | UCD405A==5)  & numberdiscussants >= 1] <- 3
discussant.1.vote[is.na(discussant.1.vote)] <- 0
#
discussant.2.vote <- NA
discussant.2.vote[UCD405B==1 & numberdiscussants >= 2] <- 1
discussant.2.vote[UCD405B==2 & numberdiscussants >= 2] <- 2
discussant.2.vote[(UCD405B==3 | UCD405B==5) & numberdiscussants >= 2] <- 3
discussant.2.vote[is.na(discussant.2.vote)] <- 0
#
discussant.3.vote <- NA
discussant.3.vote[UCD405C==1  & numberdiscussants >= 3] <- 1
discussant.3.vote[UCD405C==2  & numberdiscussants >= 3] <- 2
discussant.3.vote[(UCD405C==3 | UCD405C==5)  & numberdiscussants >= 3] <- 3
discussant.3.vote[is.na(discussant.3.vote)] <- 0
#
number.clinton.voters <- as.numeric(discussant.1.vote==1) + as.numeric(discussant.2.vote==1) + as.numeric(discussant.3.vote==1)
number.trump.voters <- as.numeric(discussant.1.vote==2) + as.numeric(discussant.2.vote==2) + as.numeric(discussant.3.vote==2)
number.other.voters <- as.numeric(discussant.1.vote==3) + as.numeric(discussant.2.vote==3) + as.numeric(discussant.3.vote==3)
#
prop.clinton.voters <- number.clinton.voters / numberdiscussants
prop.trump.voters <- number.trump.voters / numberdiscussants
prop.other.voters <- number.other.voters / numberdiscussants
#
wtd.table(presvote, prop.clinton.voters, weights=weight)
wtd.table(presvote, prop.trump.voters, weights=weight)
wtd.table(presvote, prop.other.voters, weights=weight)
#
#
# POLITICAL SOPHISTICATION
#
voted2012 <- CC16_316
voted2012[CC16_316==4] <- 1
voted2012[CC16_316!=4 | is.na(CC16_316)] <- 0
#
romneyvote2012 <- CC16_326
romneyvote2012[romneyvote2012 > 2] <- NA
#
voted2016primary <- NA
voted2016primary[CC16_327==1] <- 1
voted2016primary[CC16_327==2] <- 0
#
primary2016vote <- NA # 1=Clinton, 2=Sanders, 3=Trump, 4=Republican not Trump
primary2016vote[CC16_328==1] <- "Clinton"
primary2016vote[CC16_328==2] <- "Sanders"
primary2016vote[CC16_328==4] <- "Trump"
primary2016vote[CC16_328>=5 & CC16_328<=8] <- "Other GOP"
primary2016vote <- as.factor(primary2016vote)
#
voted2016 <- NA
voted2016[CC16_401!=5] <- 0
voted2016[CC16_401==5] <- 1
#
news_interest <- newsint
news_interest[news_interest > 4] <- NA
news_interest <- -1 * (news_interest - 5)
news_interest <- rescale01(news_interest)
#
mediause_blog <- NA
mediause_blog[CC16_300_1==1] <- 2
mediause_blog[CC16_300_1==2] <- 1
mediause_blog <- rescale01(mediause_blog)
#
mediause_newspaper <- NA
mediause_newspaper[CC16_300_3==1] <- 2
mediause_newspaper[CC16_300_3==2] <- 1
mediause_newspaper <- rescale01(mediause_newspaper)
#
mediause_socialmedia <- NA
mediause_socialmedia[CC16_300_5==1] <- 2
mediause_socialmedia[CC16_300_5==2] <- 1
mediause_socialmedia <- rescale01(mediause_socialmedia)
#
numbersocialmediapoliticalacts <- NA
numbersocialmediapoliticalacts <- 1 + ((-1 * CC16_300d_1 + 2) + (-1 * CC16_300d_2 + 2) + (-1 * CC16_300d_3 + 2) + (-1 * CC16_300d_4 + 2) + (-1 * CC16_300d_5 + 2))
numbersocialmediapoliticalacts <- rescale01(numbersocialmediapoliticalacts)
#
knowledge_correcthouse <- NA
knowledge_correcthouse[CC16_321a==1] <- 2
knowledge_correcthouse[CC16_321a!=1] <- 1
#
knowledge_correctsenate <- NA
knowledge_correctsenate[CC16_321b==1] <- 2
knowledge_correctsenate[CC16_321b!=1] <- 1
#
knowledge_correcthousesenate <- NA
knowledge_correcthousesenate[knowledge_correcthouse==2 & knowledge_correctsenate==2] <- 2
knowledge_correcthousesenate[knowledge_correcthouse==1 & knowledge_correctsenate==1] <- 1
knowledge_correcthousesenate[knowledge_correcthouse==1 & knowledge_correctsenate==2] <- 1
knowledge_correcthousesenate[knowledge_correcthouse==2 & knowledge_correctsenate==1] <- 1
knowledge_correcthousesenate <- rescale01(knowledge_correcthousesenate)
#
demplacement <- CC16_340g
demplacement[demplacement > 7] <- NA
#
gopplacement <- CC16_340h
gopplacement[gopplacement > 7] <- NA
#
knowledge.ideologyparties <- NA
knowledge.ideologyparties[CC16_340g==8 | is.na(CC16_340g)] <- 0
knowledge.ideologyparties[CC16_340h==8 | is.na(CC16_340h)] <- 0
knowledge.ideologyparties[demplacement >= gopplacement] <- 0
knowledge.ideologyparties[demplacement < gopplacement] <- 1
#
politicalsophisticationmat <- data.frame(knowhouse = knowledge_correcthouse - 1,
	knowsenate = knowledge_correctsenate - 1,
	knowideology = knowledge.ideologyparties)
#
psych::alpha(politicalsophisticationmat)
politicalsophistication.scores <- rowMeans(politicalsophisticationmat, na.rm=TRUE)
politicalsophistication.binary <- NA
politicalsophistication.binary[politicalsophistication.scores < 1] <- "Low"
politicalsophistication.binary[politicalsophistication.scores >= 1] <- "High"
politicalsophistication.binary <- factor(politicalsophistication.binary, levels=c("Low","High"))
#
#library(foreign)
#write.dta(data.frame(politicalknowledgescores=politicalsophistication.scores), file="c:/Dropbox/Networks/analysis/CCES2016knowledge.dta")
#
political.interest <- news_interest
#
networkmix <- NA
networkmix[prop.clinton.voters==1] <- "Only Clinton"
networkmix[prop.clinton.voters > 0 & prop.other.voters > 0 & prop.trump.voters==0] <- "Clinton + not sure"
networkmix[prop.clinton.voters > 0 & prop.trump.voters > 0] <- "Clinton + Trump"
networkmix[prop.clinton.voters==0 & prop.other.voters > 0 & prop.trump.voters > 0] <- "Trump + not sure"
networkmix[prop.trump.voters==1] <- "Only Trump"
#
networkmix <- factor(networkmix, levels=c("Only Clinton", "Clinton + not sure", "Clinton + Trump", "Trump + not sure", "Only Trump"))
#
#
round((table(networkmix) / sum(table(networkmix))), 2)
table(networkmix, politicalsophistication.binary)
colSums(table(networkmix, politicalsophistication.binary))
round((table(networkmix, politicalsophistication.binary) /
	colSums(table(networkmix, politicalsophistication.binary))), 2)
#
round((wtd.table(networkmix, weights=weight) / sum(wtd.table(networkmix, weights=weight))), 2)
wtd.table(networkmix, politicalsophistication.binary, weights=weight)
colSums(wtd.table(networkmix, politicalsophistication.binary, weights=weight))
round((wtd.table(networkmix, politicalsophistication.binary, weights=weight) /
	colSums(wtd.table(networkmix, politicalsophistication.binary, weights=weight))), 2)
#
#
countytype <- NA
countytype[per_gop_2016 < .4] <- "Blue"
countytype[per_gop_2016 >= .4 & per_gop_2016 <= .6] <- "Purple"
countytype[per_gop_2016 > .6] <- "Red"
#
df <- data.frame(presvote, networkmix, politicalsophistication.binary, weight, countytype)
df2 <- na.omit(df)
#
df2$primary2016vote <- factor(primary2016vote[complete.cases(df)], levels=c("Clinton", "Sanders", "Other GOP", "Trump"))
#
df3 <- ddply(df2,.(presvote),
    function(x) with(x,
      data.frame(100*wtd.table(networkmix, weights=weight)/sum(wtd.table(networkmix, weights=weight)))))
#
names(df3)[names(df3)=="Var1"] <- "Network"
names(df3)[names(df3)=="Freq"] <- "Percentage"
#
#
#pdf("Dropbox/networks/images/fig1.pdf", height=7, width=7)
#
ggplot(df3, aes(x=presvote, y=Percentage, fill=Network)) + ylim(0,100.1) +
geom_bar(stat="identity") + scale_fill_viridis(discrete=TRUE) + xlab("2016 presidential vote intention") +
ggtitle("Composition of voters' discussion networks\n2016 Cooperative Congressional Election Study (weighted)") +
theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
#
# Primary vote
#
df4 <- ddply(df2,.(primary2016vote),
    function(x) with(x,
      data.frame(100*wtd.table(networkmix, weights=weight)/sum(wtd.table(networkmix, weights=weight)))))
#
names(df4)[names(df4)=="Var1"] <- "Network"
names(df4)[names(df4)=="Freq"] <- "Percentage"
#
df4 <- na.omit(df4)
#
#pdf("Dropbox/networks/images/fig2.pdf", height=7, width=7)
#
ggplot(df4, aes(x=primary2016vote, y=Percentage, fill=Network)) + ylim(0,100.1) +
geom_bar(stat="identity") + scale_fill_viridis(discrete=TRUE) + xlab("2016 presidential primary vote") +
ggtitle("Composition of voters' discussion networks\n2016 Cooperative Congressional Election Study (weighted)") +
theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
#
# Political sophistication
#
df5 <- ddply(df2,.(presvote, politicalsophistication.binary),
    function(x) with(x,
      data.frame(100*wtd.table(networkmix, weights=weight)/sum(wtd.table(networkmix, weights=weight)))))
#
names(df5)[names(df5)=="politicalsophistication.binary"] <- "sophistication"
names(df5)[names(df5)=="Var1"] <- "Network"
names(df5)[names(df5)=="Freq"] <- "Percentage"
#
df5 <- na.omit(df5)
#df5$sophistication <- factor(df5$sophistication, levels=c("Low", "Mid", "High"))
#
#
#pdf("Dropbox/networks/images/fig3.pdf", height=7, width=7)
#
ggplot(df5, aes(x=sophistication, y=Percentage, fill=Network)) + ylim(0,100.1) + facet_wrap(~presvote) +
geom_bar(stat="identity") + scale_fill_viridis(discrete=TRUE) + xlab("Political information") +
ggtitle("Composition of voters' discussion networks\n2016 Cooperative Congressional Election Study (weighted)") +
theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
# County partisanship
#
df6 <- ddply(df2,.(presvote, countytype),
    function(x) with(x,
      data.frame(100*wtd.table(networkmix, weights=weight)/sum(wtd.table(networkmix, weights=weight)))))
#
names(df6)[names(df6)=="Var1"] <- "Network"
names(df6)[names(df6)=="Freq"] <- "Percentage"
#
df6 <- na.omit(df6)
#
#pdf("Dropbox/networks/images/fig4.pdf", height=7, width=7)
#
ggplot(df6, aes(x=countytype, y=Percentage, fill=Network)) + ylim(0,100.1) + facet_wrap(~presvote) +
geom_bar(stat="identity") + scale_fill_viridis(discrete=TRUE) + xlab("County type") +
ggtitle("Composition of voters' discussion networks\n2016 Cooperative Congressional Election Study (weighted)") +
theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
#
#  BAYESIAN ALDRICH-MCKELVEY SCALING
#
selfplacement <- CC16_340a
selfplacement[selfplacement < 1 | selfplacement > 7] <- NA
#
placements <- data.frame(CC16_340c,CC16_340d, CC16_340e, CC16_340f,
	CC16_340g, CC16_340h, CC16_340i)
placements[placements < 1 | placements > 7] <- NA
names(placements) <- c("Obama","Clinton","Trump","Garland",
	"DemParty","RepParty","SupCourt")
#
combined.placements <- data.frame(self=selfplacement, placements)
#
cutoff <- 3   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(placements)) >= cutoff & apply(placements, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- combined.placements[index,]
#
T <- T - 4
#
self <- T[,1]
TT <- T[,-1]
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
zhat[1] ~ dnorm(0,1)T(-1.1,-0.9)
zhat[2] ~ dnorm(0,1)T(,0)
zhat[3] ~ dnorm(0,1)T(0.9,1.1)
for (j in 4:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(-1,-1,1,rnorm(q-3))
)}
#
#BAM.sim <- jags.model(textConnection(modelstring),
#    data = list('z' = z, 'q' = q, 'N' = N),
#    inits = inits, n.chains = 2, n.adapt = 2500)
#
#update(BAM.sim, n.iter=2500)
#
params <- c("zhat", "a", "b")
#samps <- coda.samples(BAM.sim, params, n.iter=5000, thin=5)
#
# Write results to disk:
#save(samps,file="Dropbox/networks/analysis/CCES2016_samps.Rda")
#
# Load results from disk:
load("Dropbox/networks/analysis/CCES2016_samps.Rda")
#
libcon_a <- samps[,grep(pattern="a[", colnames(samps[[1]]), fixed=TRUE)]
libcon_b <- samps[,grep(pattern="b[", colnames(samps[[1]]), fixed=TRUE)]
libcon_all.a <- do.call("rbind", libcon_a)
libcon_all.b <- do.call("rbind", libcon_b)
#
alpha.mu <- summary(libcon_a)[[1]][,1]
beta.mu <- summary(libcon_b)[[1]][,1]
#
# Note: we need to map back to original data
#
alpha.means <- NA
alpha.means[index] <- alpha.mu
#
# Calculate ideal points
#
libcon_self <- self
#
libcon_nsamp <- nrow(libcon_all.a)
libcon_nresp <- ncol(libcon_all.a)
#
libcon_idealpt <- rep(NA, libcon_nsamp*libcon_nresp)
dim(libcon_idealpt) <- c(libcon_nsamp,libcon_nresp)
for (i in 1:libcon_nresp){
for (j in 1:libcon_nsamp){
libcon_idealpt[j,i] <- ((libcon_self[i] - libcon_all.a[j,i]) / libcon_all.b[j,i])
}}
#
libcon_idealpt.percentiles <- rep(NA,libcon_nresp*3)
dim(libcon_idealpt.percentiles) <- c(libcon_nresp,3)
colnames(libcon_idealpt.percentiles) <- c("5%","50%","95%")
for (i in 1:libcon_nresp){
libcon_idealpt.percentiles[i,] <- quantile(libcon_idealpt[,i],probs=c(0.05,0.5,0.95),na.rm=TRUE)
}
libcon_idealpt.estimates <- libcon_idealpt.percentiles[,2]
#
bamidealpoints <- NA
bamidealpoints[index] <- libcon_idealpt.estimates
#
#
#
#
df <- data.frame(alpha.means, bamidealpoints, party=factor(partyid, labels=c("D","I","R")), networkmix, weight)
df2 <- na.omit(df)
#
#dfAtt<-ddply(df2,~party+networkmix,function(x) mean_cl_boot(x$alpha.means, conf.int=.8))
#
#
A1 <- ggplot(df2, aes(x=networkmix, y=alpha.means, fill=networkmix, weight=weight)) +
geom_boxplot() +
#scale_fill_brewer(palette="Spectral") +
scale_fill_viridis(discrete=TRUE, option="viridis") +
xlab("Network Composition") +
ylab(expression(paste("Distortion Parameter (", alpha[i], ")"))) +
ggtitle("Ideological Bias\n") +
guides(fill=guide_legend(title="")) +
theme(plot.title = element_text(hjust = 0.5),
#		axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
theme(legend.position="none")
#theme(plot.title = element_text(hjust = 0.5), legend.position="none")
#
A2 <- ggplot(df2, aes(x=networkmix, y=bamidealpoints, fill=networkmix, weight=weight)) +
geom_boxplot() +
#scale_fill_brewer(palette="Spectral") +
scale_fill_viridis(discrete=TRUE, option="viridis") +
xlab("Network Composition") +
ylab(expression(paste("Corrected Liberal-Conservative Self-Placement (", x[i], ")"))) +
ylim(-5,5) +
ggtitle("Ideological Ideal Point\n") +
guides(fill=guide_legend(title="")) +
theme(plot.title = element_text(hjust = 0.5),
#		axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
#theme(plot.title = element_text(hjust = 0.5), legend.position="none")
#
#pdf("Dropbox/networks/images/BAM_alphaidealpoints_networks.pdf", height=5, width=10)
grid.arrange(A1, A2, ncol=2, widths=c(4.1,5.9))
dev.off()
#
#
#
#
# DISTORTION PARAMETER
#
temp.df <- data.frame(scale(alpha.means), relevel(networkmix, ref = 3),
	rescale01(partyid7), rescale01(ideology), rescale01(income), rescale01(age), rescale01(education),
	rescale01(female), rescale01(black), rescale01(latino), rescale01(unionmember),
	rescale01(religiosity.scores), county.gop2016, politicalsophistication.scores)
colnames(temp.df) <- c("alphameans", "networkmix",
	"partyid", "ideoselfplace", "income", "age", "education",
	"female", "black", "latino", "unionmember",
	"religiosity", "countygop2016", "politicalsophistication")
#
temp.df2 <- na.omit(temp.df)
#
#temp.df2$politicalsophistication <- factor(ntile(temp.df2$politicalsophistication, 2),
#   labels=c("Below Median","Above Median"))
#
temp.df2$politicalsophistication3[temp.df2$politicalsophistication < 0.34] <- 1
temp.df2$politicalsophistication3[temp.df2$politicalsophistication > 0.34 & temp.df2$politicalsophistication < 0.68] <- 2
temp.df2$politicalsophistication3[temp.df2$politicalsophistication > 0.68] <- 3
#
networkmix.binary <- model.matrix(~ (temp.df2$networkmix + temp.df2$politicalsophistication)^2)
temp.df2$onlyclinton <- networkmix.binary[,2]
temp.df2$clintonnotsure <- networkmix.binary[,3]
temp.df2$trumpnotsure <- networkmix.binary[,4]
temp.df2$onlytrump <- networkmix.binary[,5]
temp.df2$onlyclinton.sophistication <- networkmix.binary[,7]
temp.df2$clintonnotsure.sophistication <- networkmix.binary[,8]
temp.df2$trumpnotsure.sophistication <- networkmix.binary[,9]
temp.df2$onlytrump.sophistication <- networkmix.binary[,10]
#
temp.df2$partyidstrength <- NA
temp.df2$partyidstrength <- abs(temp.df2$partyid - 0.5)*2
#
onlyclinton.interflexplot <- inter.binning(Y = "alphameans", D = "onlyclinton", X = "politicalsophistication", Z = c("partyid", "ideoselfplace",
	"income", "age", "education", "female", "black", "latino", "unionmember", "religiosity", "countygop2016",
	"clintonnotsure","trumpnotsure","onlytrump","clintonnotsure.sophistication","trumpnotsure.sophistication","onlytrump.sophistication"),
	cutoffs=c(0,0.25,0.5,0.75,1), data = temp.df2, Ylabel = "Alpha", Dlabel = "Only Clinton", Xlabel="Political Knowledge")
#
clintonnotsure.interflexplot <- inter.binning(Y = "alphameans", D = "clintonnotsure", X = "politicalsophistication", Z = c("partyid", "ideoselfplace",
	"income", "age", "education", "female", "black", "latino", "unionmember", "religiosity", "countygop2016",
	"onlyclinton","trumpnotsure","onlytrump","onlyclinton.sophistication","trumpnotsure.sophistication","onlytrump.sophistication"),
	cutoffs=c(0,0.25,0.5,0.75,1), data = temp.df2, Ylabel = "Alpha", Dlabel = "Clinton + Not Sure", Xlabel="Political Knowledge")
#
trumpnotsure.interflexplot <- inter.binning(Y = "alphameans", D = "trumpnotsure", X = "politicalsophistication", Z = c("partyid", "ideoselfplace",
	"income", "age", "education", "female", "black", "latino", "unionmember", "religiosity", "countygop2016",
	"onlyclinton","clintonnotsure","onlytrump","onlyclinton.sophistication","clintonnotsure.sophistication","onlytrump.sophistication"),
	cutoffs=c(0,0.25,0.5,0.75,1), data = temp.df2, Ylabel = "Alpha", Dlabel = "Trump + Not Sure", Xlabel="Political Knowledge")
#
onlytrump.interflexplot <- inter.binning(Y = "alphameans", D = "onlytrump", X = "politicalsophistication", Z = c("partyid", "ideoselfplace",
	"income", "age", "education", "female", "black", "latino", "unionmember", "religiosity", "countygop2016",
	"onlyclinton","clintonnotsure","trumpnotsure","onlyclinton.sophistication","clintonnotsure.sophistication","trumpnotsure.sophistication"),
	cutoffs=c(0,0.25,0.5,0.75,1), data = temp.df2, Ylabel = "Alpha", Dlabel = "Only Trump", Xlabel="Political Knowledge")
#
A <- inter.plot(onlyclinton.interflexplot,
	ylab="Marginal Effect on alpha", xlab="Moderator: Political Knowledge", main="Network: Only D\n")
B <- inter.plot(clintonnotsure.interflexplot,
	ylab="Marginal Effect on alpha", xlab="Moderator: Political Knowledge", main="Network: D + Not Sure\n")
C <- inter.plot(trumpnotsure.interflexplot,
	ylab="Marginal Effect on alpha", xlab="Moderator: Political Knowledge", main="\nNetwork: R + Not Sure\n")
D <- inter.plot(onlytrump.interflexplot,
	ylab="Marginal Effect on alpha", xlab="Moderator: Political Knowledge", main="\nNetwork: Only R\n")
#
#pdf("Dropbox/networks/images/interflex_CCES2016.pdf", height=10, width=10)
#
ggarrange(A, B, C, D,  labels = c("","","",""), ncol = 2, nrow = 2)
#
dev.off()
#
#
alpha.lm.all <- lm(alphameans ~ networkmix +
	partyid + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2016 + politicalsophistication + politicalsophistication:networkmix,
	data=temp.df2)
#
apsrtable(alpha.lm.all)
#
#
# Interaction tests
#
alpha.lm.all.inttest <- lm(alphameans ~ networkmix +
	partyid + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2016 +
	factor(politicalsophistication3==3) + factor(politicalsophistication3==3):partyid + factor(politicalsophistication3==3):ideoselfplace + factor(politicalsophistication3==3):income + factor(politicalsophistication3==3):age + factor(politicalsophistication3==3):education +
	factor(politicalsophistication3==3):female + factor(politicalsophistication3==3):black + factor(politicalsophistication3==3):latino + factor(politicalsophistication3==3):unionmember +
	factor(politicalsophistication3==3):religiosity + factor(politicalsophistication3==3):countygop2016 + factor(politicalsophistication3==3):networkmix,
	data=temp.df2)
#
summary(alpha.lm.all.inttest)
#
#
alpha.lm.low <- lm(alphameans ~ networkmix +
	partyid + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2016,
	data=temp.df2[temp.df2$politicalsophistication3<=2,])
#
alpha.lm.high <- lm(alphameans ~ networkmix +
	partyid + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2016,
	data=temp.df2[temp.df2$politicalsophistication3==3,])
#
apsrtable(alpha.lm.low, alpha.lm.high)
#

#
# Check partisan strength
#
alpha.lm.low2 <- lm(alphameans ~ networkmix +
	partyid + partyidstrength + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2016,
	data=temp.df2[temp.df2$politicalsophistication3<=2,])
#
alpha.lm.high2 <- lm(alphameans ~ networkmix +
	partyid + partyidstrength + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2016,
	data=temp.df2[temp.df2$politicalsophistication3==3,])
#
apsrtable(alpha.lm.low2, alpha.lm.high2)
#


#
# Equality of coefficients tests
#
equalitycoefs <- function(b1,sd1,b2,sd2){
    numerator <- b1 - b2
    denominator <- sqrt(((sd1^2) + (sd2^2)))
    zscore <- numerator / denominator
    pvalue <- (2*pnorm(abs(zscore), lower.tail = F))
    out <- c(numerator, denominator, zscore, pvalue)
    names(out) <- c("numerator", "denominator", "zscore", "pvalue")
    return(out)
}
#
equalitycoefs(b1=0.04, sd1=0.10, b2=0.14, sd2=0.06)
equalitycoefs(b1=0.17, sd1=0.15, b2=0.18, sd2=0.08)
equalitycoefs(b1=-0.14, sd1=0.18, b2=-0.24, sd2=0.09)
equalitycoefs(b1=-0.18, sd1=0.10, b2=-0.38, sd2=0.06)
#
equalitycoefs(b1=-0.18, sd1=0.10, b2=0.04, sd2=0.10)
equalitycoefs(b1=-0.38, sd1=0.06, b2=0.14, sd2=0.06)
#
equalitycoefs(b1=-0.22, sd1=0.14, b2=-0.52, sd2=0.08)
#
#
# Are weights needed? No correlation = non-informative sampling (Pfeffermann & Sverchkov, 1999)
#
cor(alpha.lm.low$residuals, weight[complete.cases(temp.df)][temp.df2$politicalsophistication3<=2])
cor(alpha.lm.high$residuals, weight[complete.cases(temp.df)][temp.df2$politicalsophistication3==3])
#
apsrtable(alpha.lm.low, alpha.lm.high, stars=1, lev=c(0.1))
apsrtable(alpha.lm.low, alpha.lm.high, stars=1, lev=c(0.05))
apsrtable(alpha.lm.low, alpha.lm.high, stars=1, lev=c(0.01))
#
#
alpha.lm.low.weights <- lm(alphameans ~ networkmix +
	partyid + partyidstrength + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2016,
	data=temp.df2[temp.df2$politicalsophistication3<=2,],
	weights=weight[complete.cases(temp.df)][temp.df2$politicalsophistication3<=2])
#
alpha.lm.high.weights <- lm(alphameans ~ networkmix +
	partyid + partyidstrength + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2016,
	data=temp.df2[temp.df2$politicalsophistication3==3,],
	weights=weight[complete.cases(temp.df)][temp.df2$politicalsophistication3==3])
#
apsrtable(alpha.lm.low.weights, alpha.lm.high.weights, stars=1, lev=c(0.1))
apsrtable(alpha.lm.low.weights, alpha.lm.high.weights, stars=1, lev=c(0.05))
apsrtable(alpha.lm.low.weights, alpha.lm.high.weights, stars=1, lev=c(0.01))
#
#
X <- na.omit(data.frame(alpha.means, party=factor(partyid, labels=c("Democrat","Independent","Republican")), networkmix))
#
#plot(allEffects(combined.lm), multiline=TRUE)
#
XX <- aggregate(alpha.means ~ party + networkmix, data=X, mean)
#XX <- aggregate(alpha.means ~ party + networkmix, data=X, function(x){sd(x, na.rm=TRUE)/sqrt(length(x))})
#
#pdf("Dropbox/networks/images/BAM_alphameans_networks.pdf", height=5, width=6.5)
#
ggplot(data=XX, aes(x=party,y=alpha.means,group=networkmix,color=networkmix)) +
	geom_line(lwd=1.5) +
	scale_color_viridis(discrete=TRUE, option="viridis") +
	xlab("\nParty Identification") +
	ylab(expression(paste("Distortion Parameter (", alpha[i], ")"))) +
	ylim(-0.8,0.8) +
	ggtitle("Ideological Bias (Group Means)\n") +
	guides(color=guide_legend(title="Network Composition")) +
	theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
#
#
#
alphadf.CCES2016 <- data.frame(alpha.means, networkmix, politicalsophistication.binary, weight, partyid7)
alphadf.CCES2016$politicalsophistication.binary <- factor(alphadf.CCES2016$politicalsophistication.binary, labels = c("Low Knowledge", "High Knowledge"))
#
save(alphadf.CCES2016, file="Dropbox/networks/analysis/alphadfCCES2016.rda")
#load("Dropbox/networks/analysis/alphadfCCES2016.rda")
#
biasplot2016 <- ggplot(na.omit(alphadf.CCES2016), aes(x=networkmix, y=alpha.means, fill=networkmix, weight=weight)) +
geom_boxplot() +
facet_wrap(~politicalsophistication.binary) +
#scale_fill_brewer(palette="Spectral") +
scale_fill_viridis(discrete=TRUE, option="viridis") +
xlab("Network Composition") +
ylab(expression(paste("Distortion Parameter (", alpha[i], ")"))) +
ggtitle("Ideological Bias (CCES 2016)\n") +
guides(fill=guide_legend(title="Network Composition")) +
theme(plot.title = element_text(hjust = 0.5),
#		axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
#theme(legend.position="none")
#theme(plot.title = element_text(hjust = 0.5), legend.position="none")
#
biasplot2016
#
